﻿using System;

namespace SoftwareConsulting.BI11.NumericalMethods.DhbFunctionEvaluation
{
    /// <summary>
    /// Класс гамма-функции (интеграл Эйлера)
    /// </summary>
    public sealed class GammaFunction
    {
        private static double _sqrt2Pi = Math.Sqrt(2 * Math.PI); //корень квадратный из двух Pi
        private static double[] _coefficients = { 
                                                    76.18009172947146,
                                                    -86.50532032941677,
                                                    24.01409824083091,
                                                    -1.231739572450155,
                                                    0.1208650973866179e-2,
                                                    -0.5395239384953e-5
                                                };

        /// <summary>
        /// Возвращает значение бета-функции от указанных аргументов
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        public static double Beta(double x, double y)
        {
            return Math.Exp(LogGamma(x) + LogGamma(y) - LogGamma(x + y));
        }        

        /// <summary>
        /// Возвращает значение гамма-функции от указанного аргумента
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double Gamma(double x)
        {
            return x > 1
                        ? Math.Exp(LeadingFactor(x)) * Series(x) * _sqrt2Pi / x
                        : (x > 0 ? Gamma(x + 1) / x
                                        : double.NaN);
        }

        /// <summary>
        /// Возвращает ведущий фактор для указанног аргумента
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double LeadingFactor(double x)
        {
            double temp = x + 5.5;
            return Math.Log(temp) * (x + 0.5) - temp;
        }

        /// <summary>
        /// Возвращает значение логарифма бета-функции от указанных аргументов
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        public static double LogBeta(double x, double y)
        {
            return LogGamma(x) + LogGamma(y) - LogGamma(x + y);
        }

        /// <summary>
        /// Возвращает значение логарифма гамма-функции от указанного аргумента
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double LogGamma(double x)
        {
            return x > 1
                        ? LeadingFactor(x) + Math.Log(Series(x) * _sqrt2Pi / x)
                        : (x > 0 ? LogGamma(x + 1) - Math.Log(x)
                                        : Double.NaN);
        }
        
        /// <summary>
        /// Возвращает значение ряда в формуле Ланкзоса при указанном аргументе
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Series(double x)
        {
            double result = 1.000000000190015d;
            double term = x;
            for (int i = 0; i < 6; i++)
            {
                term += 1;
                result += _coefficients[i] / term;
            }
            return result;
        }
    }
}
